Collisional Cooling of a Charged Granular Medium 
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The dissipation rate due to inelastic collisions between 
equally charged, insulating particles in a dilute granular 
medium is calculated. It is equal to the known dissipation rate 
for uncharged granular media multiplied by a Boltzmann-like 
factor, that originates from Coulomb interactions. We include 
particle correlations by introducing an effective potential, that 
replaces the bare Coulomb potential in the Boltzmann factor. 
All results are confirmed by computer simulations. 



Submitted to PRE 



I. INTRODUCTION 

The particles in most granular materials carry a net 
electrical charge. This charge emerges naturally due to 
contact electrification during transport or is artificially 
induced in industrial processes. It is well known 
for instance, that particles always charge when trans- 
ported through a pipe. In industry, contact electrifica- 
tion is used for dry separation of different plastic materi- 
als or salts which tend to get oppositely charged and 
hence are deflected into opposite directions when falling 
through a condenser. Another application is powder var- 
nishing, where uniformly charged pigment particles are 
blown towards the object to be painted, which is oppo- 
sitely charged. 

Whereas the dynamics of electrically neutral grains 
have been studied in great detail, little is known about 
what will change, if the grains are charged. In this paper 
we present the answer for collisional cooling, a basic phe- 
nomenon, which is responsible for many of the remark- 
able properties of dilute granular media. By collisional 
cooling one means that the relative motion of the grains, 
which lets them collide and can be compared to the ther- 
mal motion of molecules in a gas, becomes weaker with 
every collision, because energy is irreversibly transferred 
to the internal degrees of freedom of the grains. 

In 1983 Haff Q showed, that the rate, at which the 
kinetic energy of the relative motion of the grains is dis- 
sipated in a homogeneous granular gas, is proportional 
to T'^/^, where T is the so called granular temperature. 
It is defined as the mean square fluctuation of the grain 
velocities divided by the space dimension: 



temperature of a freely cooling granular gas decays with 
time as t~^. We shall discuss, how these laws change, if 
the particles are uniformly charged. 

Due to the irreversible particle interactions large scale 
patterns form in granular media, such as planetary rings 
1^ or the cellular patterns in vertically vibrated granular 
layers [Q. This happens even without external driving 
|p|,p|JTl]], where one can distinguish a kinetic, a shear- 
ing and a clustering regime. The regimes depend on 
the density, the system size and on the restitution co- 
efficient Cn = —v'^/vyi, which is the ratio of the normal 
components of the relative velocities before and after a 
collision between two spherical grains. The r3/2 cool- 
ing law holds, provided the restitution coelRcient may be 
regarded as independent of Vn (ij], and the system re- 
mains approximately homogeneous The latter con- 
dition defines the kinetic regime, which is observed for the 
highest values of the restitution coefficient, whereas the 
two other regimes are more complicated because of the 
inhomogeneities. Such inhomogeneities can only occur 
as transients, if all particles are equally charged, because 
the Coulomb repulsion will homogenise the system again. 

In order to avoid additional dissipation mechanisms 
due to eddy currents within the grains we consider only 
insulating materials. Unfortunately, up to now, no con- 
sistent microscopic theory for contact electrification of 
insulators exists [ p^ . In powder processing two types 
of charge distribution are observed |^ : A bipolar charg- 
ing, where the charges of the particles in the powder can 
have opposite sign and the whole powder is almost neu- 
tral. The other case is monopolar charging, for which the 
particles tend to carry charges of the same sign and the 
countercharge is transferred to the container walls. It de- 
pends largely on the type of processing, whether one ob- 
serves bipolar or monopolar charging, which means, that 
the material of the container, the material of the powder 
and other more ambiguous things, like air humidity or 
room temperature are important ]l^ . 

The outline of this article is as follows: The next sec- 
tion specifies the model we are considering. A simple 
derivation of the dissipation rate in dilute charged gran- 
ular media based on kinetic gas theory is given in sec- 
We find, that the dissipation rate is essentially 
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A consequence of this dissipation rate is that the granular 



tion 

the one known from uncharged granular media multiplied 
with a Boltzmann factor. Section ^ compares the ana- 
lytic results with computer simulations. We find that in 
non-dilute systems the Coulomb repulsion is effectively 
reduced. This reduction will be explained, and we de- 
termine its dependence on the solid fraction of the gran- 
ular gas. In the appendix we discuss the new simula- 
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tion method we developed for this investigation. It is 
a molecular dynamics method, that avoids the so called 
brake-failure ||l4|] . 

II. THE MODEL 

In this paper, we consider monopolar charging, which 
is the usual case if insulators are transported through a 
metal pipe [^|J^. For simplicity we assume, that all par- 
ticles have the same point charge q centred in a sphere of 
diameter d and mass to. No polarisation and no charge 
transfer during contact will be considered. The parti- 
cle velocities are assumed to be much smaller than the 
velocity of light, so that relativistic effects (retardation 
and magnetic fields due to the particle motion) can be 
neglected. The electrodynamic interaction between the 
particles can then be approximated by the Coulomb po- 
tential: 

^^J=qVr,,, (2) 

where is the distance between the centers of particles 
i and j. 

We consider the collisions as being instantaneous, 
which is a good approximation for the dilute granular 
gas, where the time between collisions is much longer 
than the duration of the contact between two particles. 
As the incomplete restitution (cn < 1) is the main dis- 
sipation mechanism in granular gases. Coulomb friction 
will be neglected in this paper. Also, the dependence of 
the restitution coefficient on the relative velocity [ p^JT?! 
will be ignored, so that the constant Cn is the only mate- 
rial parameter in our model. 

The particles are confined to a volume V ~ L'^ with 
periodic boundary conditions in all three directions. The 
periodic volume can be thought of as a sufficiently homo- 
geneous subpart of a larger system, which is kept from 
expanding by reflecting walls. For vanishing particle di- 
ameter this model corresponds to the One Component 
Plasma (OCP) In the OCP a classical plasma is 

modelled by positive point charges (the ions) acting via 
the Coulomb potential, whereas the electrons are consid- 
ered to be homogeneously smeared out over the whole 
system. In the OCP the electron background cannot be 
polarised, i.e. Debye screening does not exist, as is the 
case in our model, too. 

III. ANALYTICAL RESULTS FOR DILUTE 
SYSTEMS 

In this section we derive an approximate expression for 
the dissipation rate in a dilute system of equally charged 
granular spheres, neglecting particle correlations. Basi- 
cally we apply the kinetic gas theory, but include inelastic 
collisions. Using the analytic form of the dissipation rate 



in the dilute limit derived here, we will discuss the dis- 
sipation in a non-dilute system, where correlations are 
important, in the next chapter. 

We start with calculating the collision frequency of a 
fixed particle i with any of the other particles j. If they 
were not charged, two particles would collide provided 
the relative velocity u points into the direction of the 
distance vector r = fj—fi connecting the particle centers, 
u-f > 0, and the impact parameter b — \ry. u\/u is smaller 
than the sum of the particle radii, b < 6niax — d. If the 
particles carry the charge q, they repel each other and 
the maximum impact parameter bmax becomes smaller 
than d (see Fig.|l|). By the conservation laws for angular 
momentum and for energy one gets: 

where /i ~ to/2 is the reduced mass. Eq — q^/d denotes 
the energy barrier which must be overcome to let two 
particles collide in the dilute limit. It is the difference 
of the potential energies at contact and when they are 
infinitely far apart. Eq. (|^) is independent of the actual 
form of the potential, as long as it has radial symmetry. 
(Note that energy is conserved as long as the particles do 
not touch each other.) 




FIG. 1. Particle i collides with particle j. 



Imagine a beam of particles, all having the same 
asymptotic velocity u far away from particle j. All parti- 
cles within an asymptotic cylinder of radius 5max around 
the axis through the center of j with the direction of u 
will collide with particle j. There will be irb^^^un such 
collisions per unit time, where n = N/V is the number 
density. Integrating over all relative velocities u gives the 
collision frequency of a single particle in the granular gas 
in mean field approximation: 

f^nn J d\ubl^^^{u)p{u). (4) 

U^Uq 



uo = y^2Eq/ ^ is the minimal relative velocity at infinity, 
for which a collision can occur overcoming the repulsive 
interaction. We assume that the particle velocity distri- 
bution is Gaussian with variance 3T (see (|l|)), so that 
the relative velocity will have a Gaussian distribution as 
well, with 

{u'^)=6T (5) 
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Hence, the total number of binary collisions per unit time 
and per unit volume is given by: 



Nn ^ 1/2 fn = 2^n^ (f VT ■ exp i 3. 

' niT 



(6) 



The factor 1/2 avoids double counting of collisions. This 
corresponds to textbook physics for chemical reaction 
rates as can be found for example in Present . 

Now we calculate the dissipation rate: The energy loss 
due to a single inelastic collision is: 



(7) 



where w* means the normal component of the relative 
velocity u* at the collision. It can be calculated easily 
from u*^ = u*^ — Uj^: The tangential component is de- 
termined by angular momentum conservation, 



fiub = ^u*d, 
and energy conservation gives 
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This yields 



The energy loss in one collision is therefore: 



(8) 
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Assuming a homogeneous distribution of particles, we 
eliminate the 6-dependence by averaging over the area 
TT&max (see Fig^: 
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The dissipated energy per unit time due to collisions with 
relative velocity u is then the number of such collisions 
per unit volume, 1/2 n^nb^^^ u, times the energy loss SE, 
Eq.(0). 

Finally we get the dissipation rate per unit volume in 
the dilute limit (v ~* 0) by integration over the relative 
velocity distribution: 



7=2*^^ I d-^'>^Ki,^uSE{u)p{u) 



2y/^n^d^m (1 - e^) T^/^ . exp ^- 
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The dissipation rate of an uncharged granular system in 
the dilute limit in the kinetic regime is given by H: 
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(15) 



Thus the dissipation rate ( [Til) in a monopolar charged 
dilute granular gas and the one for the uncharged case 
differ only by a Boltzmann factor, 7 = 7o-exp {—Eq/mT). 
This is the main result of the analytic treatment in this 
section. It remains valid for any repulsive pair interaction 
between the grains that has rotational symmetry. 



IV. DISSIPATION RATE FOR DENSE SYSTEMS 

In order to discuss the dissipation rate 7 in a non- 
dilute system of charged granular matter, let us recall 
the analytic form of 7 in an uncharged non-dilute sys- 
tem. The derivation is basically done by using the En- 
skog expansion of the velocity distribution function for 
dense gases . One gets the dissipation rate for a non 
dilute uncharged system: 



7 = 70 -5118(1^) 



(16) 



where 70 is given by Eq. ([l5| ) and ghs{t^) is the equi- 
librium pair distribution function of the non-dissipative 
hard-sphere fluid at contact. It only depends on the solid 
fraction v = irnd^ /6: 



5hs(i^) 



2-1/ 
2(l-i/)3 



(17) 



(Carnahan and Starling |2^, Jenkins and Richman pO|). 

Our system consists of dissipative charged hard- 
spheres (CHS). The Boltzmann factor in Eq. ( p^ is 
just the equilibrium pair distribution function at contact 
in the dilute limit for a CHS-fluid, lim^^o 5chs(i^, <?) = 
exp (—Eq/mT). So it is plausible, that the dissipation 
rate for a dense system of dissipative CHS is: 



7 = 70 • gchs{i^, q) 



(18) 



Unfortunately the literature is lacking a satisfying an- 
alytic expression for ijchs- In 1972 Palmer and Weeks |2^] 
did a mean spherical model for the CHS and derived an 
analytic expression for f/chs, but this approximation is 
poor for low densities. Many methods p3| give gchs as a 
result of integral equations, that can be solved numeri- 
cally. We do not use those approximations, but make the 
following ansatz for t^chs: 



gcY,s{v, q) ~ ghs{i^) ■ exp 



^cff(t^) 

mT 



(19) 



As in the dilute case we assume that the long range 
Coulomb repulsion modifies the pair correlation function 
of the uncharged hard sphere gas by a Boltzmann fac- 
tor. Note that the granular temperature enters the pair 
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correlation function only through this Boltzmann factor. 
The hard core repulsion is not connected with any en- 
ergy scale, so that the pair correlation function ^hs can- 
not depend on T. The effective energy barrier E^ff must 
approach Eq in the dilute limit. Hence the ansatz (19|) 
contains both the uncharged and the dilute limit, (Iq) 
respectively (p^). 

In order to check the ansatz ( |l9|) we did computer 
simulations using the MD algorithm as described in the 
appendix. Test systems of varying solid fraction v and 
particle number ranging from N — 256 to = 1024 
were prepared at a starting temperature Tq. As soon 
as the simulation starts, the granular temperature drops 
because of the inelastic collisions. We measured the dis- 
sipation rate 7 and the granular temperature during this 
evolution. According to Eq. (jl^) and Eq. jl^ ) the dissi- 
pation rate is 7 = 7o5hs(i^) ■ eTq){Ecsii')/'mT). An Ar- 
rhenius plot (ln(7/7o5hs) versus Eq/mT) should give a 
straight line whose negative slope is the effective energy 
barrier E^s- Fig. ^ shows two examples of these simu- 
lations. The Arrhenius plots are linear to a very good 
approximation. This confirms the ansatz (|l9|). Systems 
with high densities show slight deviations from linearity. 
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FIG. 2. Arrhenius-plot of the dissipation rate 7 normalised 
by the one of the uncharged system, Eq. Granular 
temperature is scaled by _Eq/m. Filled circles correspond 
to simulations of density u — 3.375 ■ 10~^ and filled squares 
u ^ 7 ■ 10"^. The linear fits yield: E^s/E^ = 0.70 for the 
lower density and Ecs/Eq — 0.27 in the other case. 

The negative slopes E^s/Eq in Fig. || are smaller 
than 1, which means, that the effective energy barrier 
is smaller than in the dilute system. The explanation 
is that two particles which are about to collide not only 
repel each other but are also pushed together by being 
repelled from all the other charged particles in the sys- 
tem. 

For dimensional reasons the effective energy barrier to 
be overcome, when two particles collide, must be of the 
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where ^ > c? is the typical distance between the charged 
particles and / is a dimensionless function. The first term 
is the Coulomb interaction Eq of the collision partners at 
contact. The second term takes the interaction with all 
other particles in the system into account. It is negative, 
because the energy barrier for the collision is reduced in 
dense systems. 

Obviously, for a dense packing, £ d, the energy bar- 
rier for a collision must vanish, i.e. Ecs\d=e = 0. More- 
over, if one takes a dense packing and reduces the radii 
of all particles infinitesimally, keeping their centers in 
place, all particles should be force free for symmetry rea- 
sons. Therefore, the energy barrier must vanish at least 
quadratically in {£ ~ d), i.e. dEcs/dd\d=i = 0. For the 
function / this implies 



and 



dfix) 



dx 



= -1. 



(21) 
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If the particle diameter d is much smaller than the typical 
distance £ between the particles, the function f{d/£) may 
be expanded to linear order. 



f{x) = Co -f cix -f 



(22) 



In linear approximation the coefficients are determined 
by ( pl| ) : Co = 2 and ci = — 1 . This determines the energy 
barrier (pO[). 

In 1969 Salpeter and Van Horn |24|j3C| ] pointed out, 
that inside a strongly coupled OCP a short-range body 
centered cubic (BCC) ordering will emerge. In the BCC 
lattice the nearest neighbour distance £ is related to the 
volume fraction v by 



d _ 2 /3 



1/3 



1.14Z.1/3. 



(23) 



Assuming a BCC structure and using the linear approx- 
imation for f{x) in (pO[), the effective energy barrier is 
therefore given by 



E^g = Eq(l~ 2.27 1/^/3 + 1.29 i/^/^) 



(24) 
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FIG. 3. The dependency of the effective energy barrier on 
the sohd fraction. Filled circles correspond to computer sim- 
ulations, the solid line is Eq. (p^. 

To test Eq. ( p^ we simulated systems with densities 
ranging from v = 0.001 to = 0.216 and determined the 
ratio Ecsiiy)/ Eq as in Fig. |[ The results are plotted in 
Fig. ^ The agreement of the theoretical formula ( |2^ ) 
with the simulations is excellent. One can sec, that in 
the dilute limit the effective energy barrier extrapolates 
to -Eq. We cannot simulate low density systems, because 
collisions are too unlikely. 

For the highest densities one cannot expect that the 
linear approximation (^2) remains valid. Also, the dense 
packing of spheres is achieved with an FCC (face cen- 
tered cubic) rather than a BCC ordering. This may be 
responsible for the systematic slight deviation from the 
theoretical curve in Fig. ||. 

The reduction of the Coulomb repulsion was also found 
in the OCP, when it was applied to dense stars [^ . 
There the analogue of the second term in ( pO| ) is called the 
"screening potential" (somewhat misleadingly, as there 
is no polarizable counter charge and hence no screening) . 
Monte Carlo simulations |^ of the OCP were interpreted 
in terms of a linear "screening potential" , which cor- 
responds to (IH) , and the analogue of the conditions ( |2l| ) 
also occurs in the plasma context , although based on 
a different physical reasoning. Corrections to the linear 
approximation are the subject of current research [ p8| . 
However, applying these more sophisticated forms of the 
"screening potential" of the OCP model to dense charged 
granular gases seems arguable as for higher densities the 
influence of the hard spheres become more and more im- 
portant and so the analogy to the OCP model, which 
uses point charges, does no longer hold. 



V. DISCUSSION 

We derived the dissipation rate of a charged granu- 
lar gas in the dilute limit. Compared to the uncharged 



case it is exponentially suppressed by a Boltzmann fac- 
tor depending on the ratio between the Coulomb barrier 
and the granular temperature. This result was obtained 
assuming a Gaussian velocity distribution, although it 
is known that in the uncharged case deviations from 
a Gaussian behaviour emerge due to the inelastic col- 
lisions | pT[ |. These deviations, however, were shown to 
have little effect on the dissipation rate |^^. As the sys- 
tem becomes less dissipative in our case, it is reasonable 
to expect that the effect of deviations from a Gaussian 
velocity distribution will be even weaker. One may say 
that a dilute granular gas with monopolar charging is 
more similar to a hard sphere gas in thermal equilibrium 
than a neutral one. 

In a dense system particle correlations enter the colli- 
sion statistics and hence the dissipation rate in two ways: 
First there is the well known Enskog correction as in the 
uncharged case. It describes that the excluded volume of 
the other particles enhances the probability that two par- 
ticles are in contact. Second the Coulomb barrier which 
colliding particles must overcome is reduced and will van- 
ish in the limit of a dense packing. 

In our simulation we did not observe the shearing or 
the clustering instability, probably because our systems 
were rather small. It is reasonable to expect, however, 
that shearing or clustering instabilities may at most exist 
as transients in the presence of monopolar charging, be- 
cause the Coulomb repulsion will homogenise the system 
in the long run. 
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APPENDIX A: 



COMPUTER SIMULATION 
METHOD 



Distinct element (or molecular dynamics (MD)) sim- 
ulations [Tsl are usually done with time step driven or 
event driven algorithms p3| . None of them is well suited 
for investigating a charged granular gas. Therefore we 
developed a new simulation scheme, which combines the 
virtues of both and will be described in this section. 

We use a "brute force" MD algorithm, which is simple 
and sufficient for our problem. More sophisticated ways 
of dealing with the long range interactions, such as the 
multipolar expansion js^], the particle-particle-particle- 
mesh 35 1 or the hypersystolic algorithms |Q should be 
used, if larger systems need to be studied. 

The event driven method for simulating the motion of 
all particles in the granular gas can be applied, whenever 
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the particle trajectories between collisions can be calcu- 
lated analytically, so that the time interval between one 
collision event and the next can be skipped in the sim- 
ulation. Obviously this is impossible in a system with 
long range Coulomb interactions. However, the idea to 
avoid the detailed resolution of a collision event in time is 
still applicable. So the velocities of the collision partners 
are simply changed instantaneously to the new values 
predicted by momentum and angular momentum conser- 
vation and an energy loss determined by the restitution 
coefficient. We shall keep this feature of event driven 
simulations. 

In the time step driven simulation method the equa- 
tions of motion of all particles in the granular gas are 
discretized using a fixed time step, which is small com- 
pared to the duration of a collision. Hence each colli- 
sion, which is modelled as an overlap between particles, 
is temporally resolved in detail. This has the advantage, 
that the formation of long lasting contacts between par- 
ticles can in principle be simulated realistically. If the 
particles carry equal charges, however, this will not hap- 
pen, so that the collisions may be approximated as be- 
ing instantaneous like in event driven simulations. Apart 
from being more efficient, this automatically avoids the 
so called brake-failure artifact (l4j , which hampers time- 
step driven molecular dynamics simulations with rapid 
relative motion. On the other hand, we need a time dis- 
cretization of the particle trajectories between collisions, 
in order to take the changing electrostatic interactions 
properly into account. 

Because of the long-range nature of the Coulomb po- 
tential, we have to include the interactions with the peri- 
odic images of the particles in the basic cell. One way to 
do this is by Ewald summation. Details of this method 
can be found in [ p^ . Another method is the minimum 
image convention: Only the nearest periodic image is 
taken into account for the calculation of the interac- 
tions. The minimum image method has the advantage, 
that it is much faster than the Ewald summation. We 
checked the validity of the minimum image method com- 
pared to the Ewald summation and found, that as long 
as Eq/mT < 10 both methods yield indistinguishable re- 
sults. This upper limit for the coupling has been found 
before in Monte-Carlo simulations of the OCP |2^. As 
our systems all satisfy this condition, we used the mini- 
mum image convention. 
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